subroutine ginzburg_landau(dt_c)
use type
use var
implicit none
integer :: i,j,k,x,label, label_l, label_r, indx
real(dp),dimension(6) :: q_tmp, q_tmp_new, dotd, delta, LdG, dQz, dQx, dQy, dQ, Qtrow,Qprow
real(dp) :: doubledot, doubledot2, nuQnu
real(dp) :: ElX, Ely, Elz, Eltot, dt_c
real(dp),dimension(:,:),allocatable :: q_new

allocate(q_new(6,tnodes))
delta(1) = 1.0d0 / 3.0d0 ; delta(4) =  1.0d0/3.0d0 ; delta(6) = 1.0d0/3.d0
delta(2) = 0.d0          ; delta(2) = 0.d0         ; delta(2) = 0.d0


do k=1,dime(3)
   do j=1,dime(2)
      do i=1,dime(1)
         label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
         if (bulk(i,j,k) == .true.) then
            q_tmp = q(:,label)
            doubledot = q(1,label)**2 + 2*q(2,label)**2 + 2*q(3,label)**2 + q(4,label)**2 + 2*q(5,label)**2 + q(6,label)**2
            call dotq (q_tmp,dotd)
            LdG = (1.d0 - UldG/3.d0)*q_tmp - ULdG * dotd + Uldg*doubledot * (q_tmp + delta )
            do x=1,6 
               
               label_l = neighbors(1,label)
               label_r = neighbors(2,label)
               Elx = ddf(1) *(q(x,label_l) - 2.d0*q_tmp(x) + q (x,label_r))
               label_l = neighbors(3,label)
               label_r = neighbors(4,label)
               Ely = ddf(2) *(q (x,label_l) - 2.d0*q_tmp(x) + q (x,label_r))
               label_l = neighbors(5,label)
               label_r = neighbors(6,label)
               Elz = ddf(3) *(q (x,label_l) - 2.d0*q_tmp(x) + q (x,label_r))

               Eltot = Elx + Ely + Elz

               q_tmp_new (x)=  q_tmp (x) - dt_c * ( LdG(x)  - Eltot )
            enddo
          elseif (boundary(i,j,k) == .true..and.infinite==.false.) then
             q_tmp = q(:,label)
             indx = lb2index(label)
             if (droplet) then
                label_l = boundary_neighbors0(1,indx)
                label_r = boundary_neighbors0(2,indx)
                dQx = (-3.0d0*q_tmp + 4.d0*q (:,label_l) - q (:,label_r)) * df2(1)
                label_l = boundary_neighbors0(3,indx)
                label_r = boundary_neighbors0(4,indx)
                dQy = (-3.0d0*q_tmp + 4.d0*q (:,label_l) - q (:,label_r)) * df2(2)
                label_l = boundary_neighbors0(5,indx)
                label_r = boundary_neighbors0(6,indx)
                dQz = (-3.0d0*q_tmp + 4.d0*q (:,label_l) - q (:,label_r)) * df2(3)
                dQ = dQx*abs(normal0(1,indx)) + dQy*abs(normal0(2,indx)) + dQz*abs(normal0(3,indx))
             elseif (channel) then
                label_l = boundary_neighbors0(5,indx)
                label_r = boundary_neighbors0(6,indx)
                dQz = (-3.0d0*q_tmp + 4.d0*q (:,label_l) - q (:,label_r)) * df2(3)
                dQ =  dQz*abs(normal0(3,indx))
             endif
             
             if (nondegenerate) then
                q_tmp_new =  q_tmp - dt_c * ( - dQ + gama(1)*(q_tmp - qb0(:,indx)))
             elseif (degenerate) then
                call surface(q_tmp,normal0(:,indx),Qtrow,Qprow,nuQnu)
                q_tmp_new =  q_tmp - dt_c * ( - dQ + 2*gama(1)*((Qtrow-Qprow)-delta*nuQnu))
             endif
           elseif (boundary_NP(i,j,k) == .true..and.NP_infinite==.false.) then
             q_tmp = q(:,label)
             indx = lb2index(label)
             label_l = boundary_neighbors1(1,indx)
             label_r = boundary_neighbors1(2,indx)
             dQx = (-3.0d0*q_tmp + 4.d0*q (:,label_l) - q (:,label_r)) * df2(1)
             label_l = boundary_neighbors1(3,indx)
             label_r = boundary_neighbors1(4,indx)
             dQy = (-3.0d0*q_tmp + 4.d0*q (:,label_l) - q (:,label_r)) * df2(2)
             label_l = boundary_neighbors1(5,indx)
             label_r = boundary_neighbors1(6,indx)
             dQz = (-3.0d0*q_tmp + 4.d0*q (:,label_l) - q (:,label_r)) * df2(3)
             dQ = dQx*abs(normal1(1,indx)) + dQy*abs(normal1(2,indx)) + dQz*abs(normal1(3,indx))
             
             if (NP_nondegenerate) then
                q_tmp_new =  q_tmp - dt_c * ( - dQ + gama(2)*(q_tmp - qb1(:,indx)))
             elseif (NP_degenerate) then
                call surface(q_tmp,normal1(:,indx),Qtrow,Qprow,nuQnu)
                q_tmp_new =  q_tmp - dt_c * ( - dQ + 2*gama(2)*((Qtrow-Qprow)-delta*nuQnu))
             endif
          else
            q_tmp_new =  q(:,label)
          endif
            q_new(:,label) = q_tmp_new
      enddo
   enddo
enddo  

q = q_new

return
contains
subroutine dotq(q,dots)
use type
implicit none
real(dp) , dimension (6):: q
real(dp),dimension(6),intent(out)  :: dots
real (dp),dimension(3,3) :: qo
integer :: i,n

qo(1,1) = q(1)
qo(1,2) = q(2)
qo(1,3) = q(3)

qo(2,1) = q(2)
qo(2,2) = q(4)
qo(2,3) = q(5)

qo(3,1) = q(3)
qo(3,2) = q(5)
qo(3,3) = q(6)

dots = 0.d0
   do i=1,3
      dots(1) =  dots(1) + qo(1,i)* qo(i,1)
      dots(2) =  dots(2) + qo(1,i)* qo(i,2)
      dots(3) =  dots(3) + qo(1,i)* qo(i,3)
      dots(4) =  dots(4) + qo(2,i)* qo(i,2)
      dots(5) =  dots(5) + qo(2,i)* qo(i,3)
      dots(6) =  dots(6) + qo(3,i)* qo(i,3)
   enddo

end subroutine dotq
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine surface(Qin,nu,Qt,Qp,nQn)
use var
use type
implicit none
real(dp), dimension(6) ,intent(in)::Qin
real(dp), dimension(3) ,intent(in)::nu
real(dp), dimension(6) ,intent(out)::Qt,Qp
real(dp), intent(out) :: nQn
real(dp), dimension(3,3) :: Qm, P, Qpp, Qtt
integer :: i,j,k,l
Qpp =0.d0;Qtt=0.d0;nQn=0.d0
Qm(1,1)=Qin(1);Qm(1,2)=Qin(2);Qm(1,3)=Qin(3);Qm(2,1)=Qin(2);Qm(2,2)=Qin(4)
Qm(2,3)=Qin(5);Qm(3,1)=Qin(3);Qm(3,2)=Qin(5);Qm(3,3)=Qin(6)

Do i=1,3
   Do j=1,3
      if (i==j) Qtt(i,j) = Qm(i,j) + Sbulk/3.d0
      if (i/=j) Qtt(i,j) = Qm(i,j) 
      
      if (i==j) P(i,j) = 1 - nu(i)*nu(j)
      if (i/=j) P(i,j) = - nu(i)*nu(j)
   enddo
enddo

Do i=1,3
   Do j=1,3
      Do k=1,3
         Do l=1,3
            Qpp(i,j) = Qpp(i,j) + P(i,k) * Qtt(k,l) * P(l,j)
         enddo
      enddo
  enddo
enddo

Do i=1,3
   Do j=1,3
      nQn = nQn + nu(i) * Qtt(i,j) * nu(j)
   enddo
enddo

Qt(1)=Qtt(1,1);Qt(2)=Qtt(1,2);Qt(3)=Qtt(1,3)
Qt(4)=Qtt(2,2);Qt(5)=Qtt(2,3);Qt(6)=Qtt(3,3)


Qp(1)=Qpp(1,1);Qp(2)=Qpp(1,2);Qp(3)=Qpp(1,3)
Qp(4)=Qpp(2,2);Qp(5)=Qpp(2,3);Qp(6)=Qpp(3,3)

endsubroutine surface


End subroutine ginzburg_landau

